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ABSTRACT 

The  representation  of  both  full  variable-step  and  fixed-leading- 
coefficient  implementations  of  BDF  formulas  using  the  modified  divided 
differences  normally  associated  with  Adams'  formulas  is  developed. 
Simple  recursion  relationships  for  the  coefficients  are  derived.   The 
technique  permits  a  simple,  compact  implementation  of  variable-step  and 
fixed-leading-coefficient  forms  of  blended  formulas. 
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J..   INTRODUCTION 

A  careful  implementation  of  modified  divided  differences  (MDDs) 
leads  to  an  efficient  implementation  of  variable-order,  variable-step 
multistep  codes  (see  Krogh  [6]  and  Shampine  and  Gordon  [8]).   These 
codes  can  be  more  efficient  than  use  of  the  Nordsieck  vector  (Gear  [2] 

and  Hindmarsh  [4] )  because  the  Nordsieck  vector  approach  uses  order  of 

2 

k  /2  additions  and  stores  in  the  predictor  phase  plus  an  additional  k 

multiplications  and  additions  and  stores  in  the  corrector  phase.   (These 
numbers  are  per  component  of  a  system.   We  will  ignore  other  operations 
that  are  independent  of  the  number  of  components.)  In  contrast,  MDDs  use 
2k  additions  and  k  multiplications  in  the  predictor  phase  followed  by  k 
additions  and  stores  in  the  corrector  phase.   To  this  must  be  added  k  - 
s  multiplications  and  stores  if  the  step  has  been  constant  for  only  s 
steps.  By  keeping  the  step  constant  whenever  possible,  MDDs  have 
considerably  less  work  per  component  for  any  but  the  smallest  values  of 
k,  and  even  if  the  step  changes  frequently  there  is  a  saving  for  modest 
to  large  k.  Although  there  is  some  additional  work  in  coefficient 
calculation  for  MDD-based  codes  making  them  perhaps  less  attractive  for 
a  very  small  number  of  equations,  those  problems  do  not  normally  use 
much  computer  time  anyway.   (Also,  if  the  step  is  constant  for  a  while, 
the  coefficients  become  constants,  saving  the  overhead.) 

One  drawback  of  MDDs  is  that  a  substantially  different  scheme  must 
be  used  for  Adams'  formulas  (for  nonstiff  equations)  to  BDF  formulas 
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(for  stiff  equations).  When  Adams'  formulas  are  used,  the  MDDs  of  y'  = 
dy/dt  are  stored,  while  when  BDF  formulas  are  used,  the  MDDs  of  y  are 
stored.   This  means  that  when  a  switch  from  a  nonstiff  to  a  stiff  method 
is  wanted  (or  vice-versa)  a  significant  disruption  is  necessary  in  the 
processing.  Nordsieck-based  codes,  on  the  otherhand,  simply  change  some 
coefficients  and  continue.  This  makes  a  "continuous  transition"  from 
nonstiff  to  stiff,  as  done  in  the  blended  methods  of  Skeel  and  Kong  [9] , 
a  reasonable  addition. 

This  paper  gives  a  technique  for  implementing  a  variable-step  form 
of  the  BDF  formulas  while  storing  the  MDDs  of  y'  (instead  of  y).   The 
terra  "variable-step  form"  is  used  in  a  technical  sense  here  to  mean  that 
the  technique  is  exactly  equivalent  to  using  a  BDF  formula  based  on 
unequal  intervals.  Formulas  based  on  unequal  intervals  have  been  shown 
to  have  generally  better  stability  properties  than  those  based  on 
constant  steps  with  interpolation  for  step  changing  (see  Gear  and  Tu 
[3]).   The  latter  method  is  the  simplest  to  use  with  the  Nordsieck 
scheme  although  the  former  have  been  used  (Byrne  and  Hindmarsh  [1]). 
The  former  method  is  more  suited  to  MDD-based  codes . 

Although  this  technique  is  more  expensive  than  storing  the  MDDs  of 
y  (by  k  multiplications  and  additions),  it  allows  a  simple  transition 
from  Adams  to  BDF  formulas  and  thus  permits  a  simple  implementation  of 
the  blended  methods.   The  additional  cost  is  not  too  significant  in  BDF 
methods  because  the  corrector  normally  requires  the  solution  of  implicit 
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equations,  and  this  is  the  most  time-consuming  part  of  the  process. 

The  paper  also  describes  the  very  simple  change  necessary  to  get  a 
fixed-leading-coefficient  form  of  the  method,  and  the  application  of 
this  form  to  a  "Blended"  version. 

The  second  section  of  this  paper  describes  the  basis  of  the 
technique.   This  is  illustrated  by  a  simple  fixed-step  second  order  BDF 
method.   The  third  section  discusses  the  calculation  of  the  BDF 
coefficients  in  a  variable-step  method  and  develops  difference  equations 
which  are  simple  to  evaluate.  The  final  section  discusses  fixed- 
leading-coefficients  and  blended  methods. 
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2.      POLYNOMIAL  INTERPOLATION  AND  UPDATING 

This  section  first  describes  the  basic  MDD  approach  to  an  Adams' 
formula-based  code  in  general  terms.   (Details  can  be  found  in  [6]  and 
[8].)  From  this  viewpoint  it  is  then  clear  how  a  BDF  formula-based  code 
can  be  implemented  storing  the  same  quantities  as  the  Adams'  scheme. 
Finally  this  Is  illustrated  with  second-order  examples. 

Suppose  that  the  solution  has  already  been  calculated  at  a  set  of  k 

points  t  ,  t   ......  t  ,,,  and  we  have  kept  the  values 

r       n   n-1       n-k+1  r 

y   y',  y'  ,,...,  y'  ,  , . .   In  the  MDD-based  method  we  actually  store  the 
J n,  Jn     yn-l      ^n-k+1  J 

1* 
modified  divided  differences  VJ  y' ,  j  *»  0, . . . ,  k-1,  which,  for  constant 

steps,  are  equal  to  the  backward  differences  of  order  j.   The  first 

process  in  calculating  the  value  at  t  , ,  ■  t  +  h  is  prediction,  and  can 
r  n+1    n 

be  written  as 


k-1 
J-0 

(1) 


'n+1    *   7n  y   n 


k_1       4* 
P^i  =  y  +  h  V   Y4   VJ   y' 
rn+l   'n     AL  n     In  n  yn 

j-0  J 


(2) 


where  the  y^   are  the  Adams-Bashforth  coefficients  (which  are  functions 
jn 

of  the  stepsize  ratios).   It  is  important  to  note  that  prediction 

consists  of  nothing  more  than  calculating  the  value  and  derivative  at 

t  .  of  the  k-th  degree  polynomial  which  matches  y  at  t  and  y'   at 

t   .,  0  <  I  <  k-1.   The  second  process  is  to  calculate  a  correction 
n-i  r 
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term.   This  depends  on  which  type  of  method,  namely  PE,  PEC,  PECE,..., 

or  solving  the  corrector  by  Iteration,  is  used.   Whatever  method  is 

used,  we  are  computing  the  value  y^,  which  is  such  that  ?'+1  -  f(y_.i) 

for  the  differential  equation  y'  »  f(y),  where  ?',,  is  defined  as  the 

derivative  of  the  polynomial  passing  through 

y  .i»  y   y' t  y'  1 1  •  • • »  y'   •   The  corrector  order  is  R+l  and  normally 
yn+l  •'n,  Jn     ^n-1      J     _. 

n-K 

R  »  k.   If  the  corrector  is  iterated  to  convergence,  $',-.    ■  f(y^_i)  anc* 

y',,  =  y'.,  •   For  the  other  methods  we  have 
yn+l   yn+l 

PE!      K+i  -  "*i  s°  1*1  -  *Vl 

1*1  •  ^.H-l' 

PEC!    Kn '  1*1  -  £(>Vi> 

PECE:    ra+l  ■  Up^) 

Finally  we  have  to  update  the  saved  Information  so  that  we  have  the 

representation  of  y  .,,  y' , , ,  y' , . . . ,  y'  ,  .  „  .  for  the  next  step  if  the 

n+1   n+1  Jn  yn-k+2-j 

order  is  the  same  (j  =  0),  increased  by  one  (j  =  1)  or  decreased  by  - j . 

The  values  of  y  .  ,  and  y'+1  are  defined  above.   The  remaining  MDDs  can 

be  updated  by  observing  that  y'  .  -  p'  .  is  the  difference  between  y'  . 

and  the  derivative  at  t   .  of  the  k-th  degree  polynomial  whose 

derivatives  are  y'    at  t   . ,  0  <  I  <  k.   Hence  this  is  exactly  the  k-th 
n-i     n-i  J 

modified  divided  difference  of  y'   ,  -1  <  1  <  k,  V  , ,  y' .   From  this  the 

•'n-j'      J     '   n+1  J 

iteration 

1      i*    1+1 
V-J    »  VJ   +  VJ^J    j  -  k-1,...,  0 
n+1    n     n+1    j      »    » 
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can  be  used,  where  the  MDD  rM  is  the  1-th  order  MDD  at  t  . ,  based  on 

n+1  n+1 

h  ■  t   .  -  t  .   It  now  has  to  be  multiplied  by  a  step-ratio  dependent 

1* 
scale  to  get  VJ   ,  the  MDD  at  t  , ,  based  on  the  next  stepslze 
n+1  n+1  r 

t  _  -  t   ..   Details  are  given  In  Krogh  or  Shamplne  and  Gordon  op_.  cit. 

After  the  update  in  Eq.  (3),  the  MDD  table     contains  one  additional 

entry,  V     and  represents  a  polynomial  through  y'   ,  -1  <  j  <  k.  This 

is  appropriate  if  the  order  is  to  be  increased  by  one.   If  not,  we  must 

discard  y'  ,  . ,  y'  ,  ., ■ . •   as  necessary — that  is  to  say  we  must  update 

the  MDD  table  so  it  represents  a  polynomial  of  minimum  degree  passing 

through  the  saved  values.  For  Adams'  method  the  saved  values  are  the 

most  recent  derivatives  so  all  that  Is  necessary  is  to  discard  the 

highest  order  MDDs  of  y'  . . 

n+1 

Now  let  us  examine  the  BDF  formulas  using  the  same  representation. 

Suppose  that  we  have  already  computed  y,y  .,...,  y  ,,  and  are 

representing  this  information  by  saving  y  and  the  MDDs  of  the 

derivatives  at  t   ,  0  <  j  <  k,  of  the  k-th  degree  polynomial  passing 

through  y   ,  0  <  1  <  k.  The  prediction  process  again  consists  of 

calculating  the  value  and  derivative  at  t   .  of  the  polynomial.   The 

predicted  values  are  given  by  the  same  formulas  (1)  and  (2)  as  the 

Adams'  code.   The  correction  process  is  also  identical  (usually 

consisting  of  solving  the  corrector  equation  starting  with  the  predicted 

value  and  using  a  quasi-Newton  iteration).  However,  note  that  we  now 

get  the  value  of  y  ,,  such  that  ?',,,  the  derivative  at  t  (1  of  the 
yn+l  ■'n+l  n+1 

polynomial  passing  through  y   . ,  -1  <   1  <  k,  approximately  satisfies 
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f    ,  ■  f(y  .,).  The  main  difference  comes  in  the  updating  process.  Now 
n+ln+1 

we  want  to  calculate  the  MDDs  V    y'   which  are  the  MDDs  of  the 

n+1  n+1 

derivatives  of  the  polynomial  passing  through  y   ,  -1  <    i  <  k-l+j  where 
j  is  +1  if  the  order  is  increased  by  one,  0  if  the  order  is  unchanged, 
and  negative  if  the  order  is  to  be  reduced. 


How  can  we  do  this?  Let  us  first  consider  the  case  of  increased 
order.   At  the  start  of  the  step  the  saved  values  represent  the  unique 

k-th  degree  polynomial  passing  through  y   ,  0  <  i  <  k.  For  such  a 

k* 
polynomial,  V   y  =  0.   At  the  end  of  the  step  the  saved  values  must 

represent  the  unique  (k+l)-st  degree  polynomial  passing  through 

y   ,  -1  <  i  <  k.  Hence,  the  polynomial  will  have  been  modified  by  the 

addition  of  a  (k+l)-st  degree  polynomial  which  is  zero  at 

t   . ,  0  <  i  <  k  and  non-zero  at  t  . ,  (for  example, 
n-i       ,  n+1         r 

k 

P  ,(t)»  IT   (t-t   )  or  any  scalar  multiple  of  it).  Hence,  the  new 
n,k      1=Q       n-l 

MDDs  at  t  are 
n 


(4) 


VJ   0'  *  VJ   +  aVJ   P'  , 
n  yn+l    n      n   n,k 

1* 

for  some  scalar  a,  where  by  VJ  9',,    we  mean  the  MDDs  at  t   of  the 

n  yn+l  n 

derivative  of  the  polynomial  through  y  ., >  y  >  y   1»»««»  y  t«   The 

coefficient  a  must  be  chosen  so  that 


Pn+1  +  aPn,k(tn+l)  ~  ^n+1  (5) 

and 
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»;+i +  ^.k'^i' a  £(w 

Finally,  we  confute  V^_  ?'    from  V^   ?'  .  using  Eq.  (3). 

This  gives  the  MDDs  when  the  order  is  raised  by  one.  What  happens 
if  it  is  to  stay  the  same  or  is  lowered?  We  will  analyze  this  by 
considering  how  to  solve  the  following  problem:   Given  the  MDDs  of  the 
derivative  of  the  k-th  degree  polynomial  through  y   ,  0  <  i  <  k,  how 
can  we  find  the  MDDs  of  the  derivative  of  the  (k-l)-th  degree  polynomial 
through  y   ,  0  <  i  <  k?  If  we  can  do  this,  we  can  apply  the  process 
repeatedly  to  lower  the  order  as  many  times  as  needed. 

The  modification  needed  consists  of  the  addition  of  a  k-th  degree 

polynomial  which  is  zero  att   .,  0  <  i  <  k  to  the  original  polynomial 

n— l 

to  reduce  its  degree  to  k-1.   Thus  we  add  a  multiple  of 


k-1 

P  ,  .  -  n   (t  -  t   .). 
n,k-l   i=Q       n-i 


The  modification  to  the  MDDs  are 


1*       -t*        i* 
VJ   y'  <-  VJ   y'  +  aVJ   P'  ,  , 
nn    nn     n   n,k-l 


(7) 


k-1* 
such  that  the  new  value  of  V     y'  is  zero.   Hence, 

n    Jn 

yJ*  y-  +  vj*  y'  -  Vj*  P'     -—2- 2—  ,   j  -  0,...,  k-1 

n   n    n  yn    n   n,k-l  yk-1*  _,         J 

n,k-l  (8) 

Thus,  the  integration  process  for  the  BDF  formula  is  as  given  below.   In 
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this  process  we  will  write  8    for  V   P'   ,  <J>  ,  for  /   6  ,  ,  and 

n,k  n       n,k        n,k  .     .      n,k 

-1  1  ^ 

<>     .    for  P     .  (t    , .)    (<+>J    ,    is   the   i-th  order  MDD   of  P'   ,    at   t    , ,    based  on 
n,  k  n,k      n+l  n,R  n,K  n+l 

h  -  t  .  .  -  t  ). 
n+l    n 

Single  Step 

(a)  Predict  using  Eqs.  (1)  and  (2). 

(b)  Correct  by  "solving" 


p'        +  ot<}>  «  f  (p  +  a<}>        ) 

n+l  n,k  rn+l  n,k 


for  a  . 
(c)  Update  MDDs  by 


(9) 


i*         i*        i 
v   f\  ,  -  V   y'  +  a6x  ,    0  <  i  <  k 
n   n+l    n   n     n,k 


k*         k 
n   n+l     n,k 


(d)  Advance  MDDs  to  t  . ,  by 

n+l 

Vk        v'        =»  Vk*  9' 
n+l   yn+l  n     yn+l 

i  i*         i+1 

v\i  y  j.i  ■  v     9'j.i  +  v  It  y\i  »  i  -  k-i,...,  o 

n+l  yn+l    n  yn+l    n+l  yn+l 

i* 

(e)  Modify  the  divided  differences  to  get  V    y'+i • 

(f)  Lower  order  if  necessary  by 


-    „k*   ,    /nk 
n    n+l  7  n+l  n+l,k 

i*  i*  i 

n+l  yn+l    n+l  yn+l    n  n+l,k 


This  order  lowering  can  be  repeated  for  k-1,  k-2,...,  as 

required. 


(10) 


(ID 


(12) 
(13) 


-  11  - 


In  fact,  (c)  and  the  last  application  of  (f)  from  the  previous  step  can 
be  combined  to  reduce  arithmetic.   This  is  described  in  the  next  section 
on  coefficient  calculation. 

Example.   2nd  order  BDF  method  with  constant  stepsize. 


and 


This  starts  with  the  values  y  , ,  y  which  are  represented  by  y 

'n-l'n  r  J    Jn 

y'    =    (y     -  y      ,)/h.      P     ,    and  its  MDDs   are  given  by 
J  n  J  n       J  n-l  n,l  G  J 

P         -   (t   -  t   )  (t   -  t     +  h) 

n,  i  n  n 

e°      -  h  ,      8*      -  a  , 

n,l  n, l 

4>      ,    ■   3h    ,      <{>      ,    =  2h 
n,l  n,l 

The  predictor  is 


P^i  -  y  +  hy'   (  =  2y  -  y  . ) 
rn+l   ;n    Jn  J  n   yn-l 


p0+i  -  "„ 


The  corrector  solves 


p;+[  +  3ho  =  «p         +  2h2a) 


or 


4  1  2h  mt         , 

yn+l   '3yn"  IVl  +  T  f  (W 


where 
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yn+l   -  Pn+1            yn+l   "   2yn  +  yn-l 
0     . - 

2ti  21i 


This  is  precisely  the  2nd  order  BDF  method.  At  the  increased  order  we 
update  the  MDDs  to  get 


followed  by 


n  yn+l   yn         2h 


_  yn-t-l  "  yn-l 
2h 

and 


vl*      a  yn+l  '  2yn  +  yn-l 
n  yn+l  h 


vl    ,     yn+l  ~  2yn  +  yn-l 
n+1  yn+l  '         h 


3yn+l  "  4yn  +  yn-l 


7°   y 
n+1  yn+l  2h 

i* 
At  fixed  stepsizes  there  is  no  further  modification  to  get  V   . ,  and  we 

can  lower  the  order  by  computing 

,0*        ,  „0  nl  ,        Q0  /Q1 


7 


n+1  yn+l         n+1   yn+l  n+1   yn+l     n+l,l'   n+1,1 


yn+l   '  yn 
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and 

Note  that  the  BDF  corrector  has  one  higher  order  than  the  predictor-. 
This  is  different  to  the  BDF  formulas  implemented  in  DIFSUB  [2]  where 
the  predictor  and  corrector  have  the  same  order. 
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3.   COEFFICIENT  CALCULATION 


Let  P  .  (t)  -  IT   (t  -  t  ).  Define  Q  ,  (t)  -  P'   (t).   It  is  easy 

n,  k      in  n  n,  k       n,  k. 

to  see  that 


«n.k(t)  "  (t  "  'n-k'  "n.k-l(t)  +  Pn.k-l(t)  (M) 

Let  the  i-th  order  divided  difference  of  Q  ,  at  t  be  w  '  .  We  have 

n,k     m     n,k 

^l  -  Q  t(t  ) 

n,k   ^n,k  m  ^^ 

M»»J+1  -  [J*'}  •\rm}'t]/(t     -  t   .  ,). 
n,k       n,k    n,k      m    m-l-i  ..  ,s 

Ultimately  we  want  to  compute  the  modified  divided  differences  of  Q  at 

t  .   First  we  will  examine  the  regular  divided  differences  w  '  .  We 
n  °  n,k 

would  like  to  generate  these  simply  so  we  look  for  a  recursion  on  k 

using  Eq.  (14).  Noting  that  P  ,  , (t  )  -  0  for  n-k+1  <  m  <  n  we  have 

n , k— 1  m 

""Ik  -  <e.  -  Vk'  Vk-i'V 

for  the  same  range  of  m.  Now  we  claim  that 


wm»j .  (t  - 1  .)  *"•£ .  +  w^*-1 

n,k     m    n-k   n,k-l    n,k-l 

-  (t   .  -  t  .  )wm»*  .  +  W*'*"1 
m-i    n-k  n,k-l    n,k-l 


(17) 


if  n-k+i  <  m  <  n  and  1  <  i. 
Proof 

For  1  ■  1  we  have  by  definition 
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W* 


,1        (tm  -  tn-k)    Qn.k-l(tm)   "   (Vl   "  Vk)    Sufe^fcl? 


n,k  t     -  t      . 

m         m-1 


Um  "     n-k'  t     -  t      , 

m         m-1 


+  Vw'Vi' 


or 


,.        ,   qn.k-l(t„>   "  ^.k-l'Vl* 

Um-1   "     n-k'  t      -  t      . 

m  m-1 


+  Vk-lV 
In  either  case   the   result   follows   from  Eq.    (15). 

Now  consider  1  >   1  and  use  induction.      Substituting  Eq.    (17)    into 


the  definition  of  w   *      we  get 

n,k         ° 


4        W*'*"1   -  W*-?'1-1 

«*»  i  ,     n»k H*J£ 

n,k  t     -  t      . 

ra  m-i 


ct  - 1  . )  *■»£-}  -  (t  ,  - 1  .)  w^?^-1 

m         n-k       n,k-l  m-i  n-k       n.k-1 


m         m-i 


Q.E.D, 
from  which  the  result  follows  directly. 


The  modified  divided  differences  are  given  by 
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n,k      nri-1    m      nrfl    m+l-i   n,k 
For  our  purposes  it  will  be  convenient  to  introduce  a  scale  of  P  ,  and 

Dy  K 

its  divided  differences  so  that  Q  ,  (t  )  »  P'  ,  (t  )  -  1.   This  is 

Ti,k  m     n,k  n 

achieved  by  dividing  by  IT  (t  -  t   .).  Hence,  we  define 

1-1  n    Brl 

0m,i  _  (tmfl  "*  ^ *  "  (tari-l  "  tmfl-i)  „m,i 

n,k  "    (t  -  t   )...(t  -  t   )     n,k 
n    n— I      n    n— k 

The  6  coefficients  defined  earlier  are  then  given  by  9    ■  0  '  .  Eq. 

LI  y   K         U.  j  K 

(17)  now  yields  the  recursion 

Qi   _  (tn-I  "  fcn-k}  6n.k-l  +  (tn+l  "  tn+l-i)  6nTk-l 

n,k  "  t  -  t  , 

n    n-k 


(18) 


1  <    i<k 
Because  of  the  scaling,  6  ,  -  1.   We  will  use  the  notation  of  Krogh 
(op»  cit«),  namely, 

Vn+1)  "  Sh-i  -  Vi 

60(n+l)  -  1 

P±(n+1)  -  B^^n+1)  C1.1(«H-l)/51_1(n) 


from  which  we   find   that 

"n.k  "  Vn+1>   "S.'k 

It  is  easy  to  show  that  W11',  -  k+1.   (It  is  the  k-th  divided  difference 

n,k 

of  the  derivative  of  a  monic  (k+l)-st  degree  polynomial.)  Hence,  we  know 

Ok  i 

9      ,    and  9        ,   k  -  0,    1 so  we   can  calculate  9  using   the 

n,  k  n,  K  n,  Ic 
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recursion  Eq.  (18)  rewritten  as 

<     IC.(n+l)  -  C.(n+1)]  6*     +  ?   (n+1)  6*~) 
e1   , J 1 n.1-1    i-1 n.1-1 

a,  J  Cj(n+1)  -  C0(n+1)' 

(19) 

for  i  -  1,...,  j-1;  j  -  2,...,  k  .   Note  that  if  the  steps izes  have  been 

constant  for  s  steps  (that  is, 

'nfl  "  fcn  '  fcn  -  Cn-1  —  -  ^2-8  "  tnfl-s)»  then  6j(n+1)  "   X  and 

e    ■  (j+D/(j+l-i)  for  j  <  s.   Thus  we  can  start  with  j  -  s  -  1  and 
n»  J 

calculate  6    explicitly  for  i  -  1,...,  j.   If  s  >  k,  the  previously 

calculated  values  of  9   ,  ,  can  be  used.   We  also  need 

n-l,k 

*n,k  '  jQ  9n,k  "  J-0.....k 

(20) 

and 


♦;!k  -  v*+i>  ek(n+i) 


(21) 


in  the  corrector  step 


Finally,  note  that  the  action  represented  by  Eq.  (13)  can  be 

delayed  until  the  start  of  the  next  step  when  Eqs.  (1),  (2),  and  (10) 

1* 
must  be  executed.   Suppose  that  VJ   y'  are  the  MDDs  left  by  Eq.  (11)  in 

the  previous  step  update  to  the  new  stepslze.   Eqs.  (1)  and  (2)  can  be 

written  as 
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*n+l   /;n   n  'n    n-1  n,k 


J-o 


J*  -'  _  a    ,.° 


(22) 


*-rt  n  ;n    n-1   n,k    n,k 


J-o 


k_1      i*  1 

p^i  -  y  +  h  I     Y4  (VJ  y'  -  6  .  9J  ) 
rn+l   'n     .^0  jn  n  'n    n-1  n,k 


k-1  k-1 

»  y  +  h  [  y.      VJ   y'  -  h6   ,   £  y4      0J 

n    j-o  Jn  n  n    n"1  j-o  Jn  n,k 

The  order  lowering  represented  by  Eq.  (13)  can  now  be  handled  by 
changing  the  first  member  of  Eq.  (10)  to 

V1*  ?'  .  -  V1*  y'  +  (a  -  6   , )  61 
n  7n+l    n  •'n        n-1   n,k 


(23) 


(24) 


Alternatively,  the  6    terms  can  be  omitted  from  Eqs.  (22)  and  (23)  if 
the  summations  are  taken  from  0  to  k.   This  corresponds  to  using  a 
predictor  of  one  order  higher — which  is  possible  if  the  order  did  not 
increase  in  the  previous  step. 

Examination  of  the  equations  given  above  reveals  that  several 
divisions  can  be  saved  by  introducing  another  normalization  of  0  so  that 
9    ■  l.   Our  program  computes  0  using  Eq.  (18)  and  then  makes  this 

11)  Iv 

final  normalization  before  computing  $  and  <J> 
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4..   FIXED  LEADING  COEFFICIENT  FORMULAS  AND  BLENDED  METHOD 

The  corrector  equation  which  must  be  solved  at  each  step  for  stiff 
equations  is 

p'  +  a<J>   «  f  (p  +  a<j>   ) 

(25) 

where  we  have  dropped  the  numerous  subscripts.  For  stiff  equations  f 

is  large  so  a  quasi-Newton  iteration  is  often  used  to  solve  Eq.  (25). 

This  means  that  the  matrix  [I  -  J<(>  /$   ]   must  be  decomposed,  where 

J  «  f  .  Unfortunately,  in  variable-step  methods  <j>   /<J>   is  only  constant 

if  the  step  and  order  have  been  fixed  for  at  least  k  steps.   This  is  in 

contrast  to  interpolatory  methods,  usually  implemented  with  Nordsieck 

vectors,  in  which  <j)   /<}>   is  a  function  of  the  current  stepsize  and  order 

only.   In  the  latter  methods,  it  is  not  necessary  to  redecompose 

[I  -  J<J>   /<J>  ]  unless  the  stepsize  or  order  changes.   In  [5],  Jackson  and 

Sacks-Davis  discuss  fixed-leading-coefficient  methods.   These  methods 

have  the  advantage,  along  with  interpolatory  methods,  that  the  ratio 

<t>  /<)>  does  not  change  unless  the  stepsize  or  order  changes.  There  are 

two  equivalent  ways  of  viewing  these  methods.   In  a  variable-step 

framework,  we  want  to  use  the  k-th  order  formula 


k+1 

n    °   n   i»i   n*1  a-1 

(25) 

where  P.  is  the  coefficient  for  the  fixed-step,  k-th  order,  k-step  BDF 

method.   This  determines  the  a    uniquely  unless  all  stepsizes  are 
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equal,  in  which  case  we  get  uniqueness  (and  continuity)  by  choosing 

a  ,  . ,  ■  0.   Another  view  of  the  method  is  to  start  with  the  unique  k-th 
n,k+l 

degree  predictor  polynomial  through  y   ,  1  <  i  <  k+1,  interpolate 

values  of  y"   .  to  a  set  of  k  equally  spaced  points  at 

t  .  -  (i-l)h,  1  <  i  <  k,  and  use  the  fixed-step,  k-th  order  BDF  formula 

k 

y   -  Bft  hy'  +  £   a.  f      . 
yn    0  ;n    L.      i  'n-i 

(26) 
It  is  simple  to  show  that  these  are  equivalent  by  considering  the 
expansion  of  J  into  linear  sums  of  the  y    to  convert  from  Eq.  (26) 

to  Eq.  (25).   This  is  a  unique  mapping  and  since  Eqs.  (25)  and  (26)  are 
uniquely  determined  by  the  k-th  order  criterion;  the  results  must  be 
identical. 


This  fixed-leading-coefficient  form  can  be  implemented  using  the 
MDD  scheme  as  follows: 

a.  Predict  using  Eqs.  (1)  and  (2)  with  sums  to  k. 

b.  Solve  the  corrector  equation 

y'  -  p'  +  a$°  -  f(p  +  a*"1)  -  f(y) 
where  4>   and  <J>   are  the  fixed-step  coefficients, 
c   Update  using  Eqs.  (24)  and  (11). 

Blended  methods  [9]  can  be  implemented  by  noting  the  difference 
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between  Adams  and  BDF  methods  arise  only  in  the  definition  of  8  and  <J>. 

B      B 
Let  us  write  9  and  <J>   for  the  coefficients  introduced  above.  The 

equivalent  Adams  coefficients  are 

eA,i  -  6<u  »    0  <  i  <  k 
n,k    ik 

where  6.,  is  the  Kronecker  delta, 
ik 

n,k 

♦  i   -  hy  . 
n,k      n,k 

k* 
where  y     ,  is  the  Adams-Bashforth  coefficient  of  7   y' .   (As  before,  we 
n,k  n  ■'n 

can  have  a  fixed-leading  coefficient  form  by  choosing  y     ,  ■  y.»  the 
fixed-step-coefficient.)  In  blended  methods  we  form  the  linear 
combination 

A  -  rJB  =»  0 

where  A  ■  0  represents  the  Adams  method,  B  -  0  represents  the  BDF 
method,  J  is  the  (approximate)  Jacobian^  and  r  is  a  real  scalar.   For 
constant  J  we  can  save  the  MDD  representation  of  [I  -  rJ]    (P.  -  rJPR) 
where  P  and  Pn  are  the  polynomials  through  the  saved  values  for  Adams 
and  BDF  methods,  respectively.  To  accomplish  this,  we  simply  solve 

p'  +  <!>A'0  a  -  r<>B'0ja  -  f(p  +  .J,*'"1  a  -  r*8'"1!  a) 

and  then  update  the  MDD's  using  Eq.  (10)  modified  to 
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V1*  y'    .  v1*  y'  +  a6A'J  -  rJa9B»J 
n  'n+1    n  Jn  n,k       n,k 

As  before,  this  can  be  modified  to  combine  steps  to  get  the  replacement 

for  Eq.  (24)  by  replacing  a  by  (a  -  6    )  In  Eq.  (27).   With  the 

It 

normalization  of  6   ,  ■  1  (for  Adams  and  BDF),  6   ,  is  precisely  the  a 

n,k  n-1    r  J 

from  the  previous  step. 
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